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We study the behavior under perturbations of different versions of Bak-Sneppen 
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(BS) model in 1+1 dimension. We focus our attention on the damage-spreading 
features of the BS model with both random as well as deterministic updating, 
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with sequential as well as parallel updating. In addition, we compute analytically 
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the asymptotic plateau reached by the distance after the growing phase. 
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I. INTRODUCTION 

A great deal of evidence has been put forward in recent years for the appearance of power 
law statistics in nature: A wide variety of phenomena, from earthquakes [|T],3 to biological 
evolution [R|-f7[ , from surface growth [[| to fluid displacement in porous media M , exhibit scale 
invariance in both space and time. To explain these all-pervading power-law tails, Bak, Tang, 
and Wiesenfeld introduced the concept of self-organized criticality (SOC) fllQf . In a nutshell, 
SOC means that certain driven spatially extended systems evolve spontaneously towards a crit- 
ical globally stationary dynamical state with no characteristic time or length scales pl| . This 
scale invariance implies that the correlation length in these systems is infinite and consequently 
a small (local) perturbation can produce a global (maybe even drastic) effect. This possibil- 
ity leads naturally to the study of the sensitivity to perturbations in (self-organized) critical 
systems. 

To study the propagation of local perturbations (damage spreading) one can borrow a tech- 
nique from dynamical systems theory. Let us consider, for instance, two copies of the same 
dynamical system (let us say, for instance, a ID map) , with slightly different initial conditions. 
By following the dynamics of both copies and studying the evolution in time of the "distance" 
d(t) between them, it is possible to quantify the effect of the initial perturbation. Indeed, 
assuming that the distance d(t) grows exponentially, and defining the Lyapunov exponent A 
via 

D(t) = D exp{\t), (1) 

three different behaviors can be distinguished, corresponding to A being either positive, negative 
or zero. The case A > corresponds to the so-called chaotic systems, where the extremely high 
sensibility to initial conditions leads to exponentially diverging trajectories on a strange (or 
chaotic) attractor. The case A < 0, instead, characterizes those systems in which the dynamics 
has an attractor and any initial perturbation is "washed out" with exponential rapidity. 
The boundary case, A = 0, admits, in turn, a whole class of functions D(t), namely 

D[t) ~ t a . (2) 

where a is some exponent, characteristic of the system. In particular, a > corresponds to 



weak sensitivity to initial conditions while a < corresponds to weak insensitivity to initial 
conditions (as an example, the reader is referred to Refs. ||12|| , where this analysis is performed 
for the logistic map at its critical point |H| ) . 



Recently fl4^|lT|, a similar analysis was performed on several versions of the Bak-Sneppen 
(BS) model ||. Originally proposed to describe ecological evolution, this model has been paid a 
great deal of attention due to its simplicity and the fact that it exhibits self-organized criticality. 

Schematically, the BS model is defined on a lattice where at each time-step one site is chosen, 
namely the one that fulfills a global constraint (minimum in some phase space). This site is 
defined as the active site. This dynamics leads to a non-Markovian process where the activity, 
i.e. the position on the lattice of the active site, jumps on the lattice following a (correlated) 
Levy Walk. Its critical properties allow us to describe its behavior under perturbations via Eq. 
(|), with 

a = 0.32. (3) 

In this paper we analyze in detail the behavior of different versions of the BS model. As 
we shall see, the particular stationary distribution of the variable does not play any role in the 
determination of the exponent a in Eq. (0). What matters is, instead, the kind of Levy Walk 
involved and the strength of the correlations. 

The paper is organized as follows. After a general description of our formalism in Sec. O, 
the BS model is discussed in Sec. |ITT[ In Sec. [TV] we analyze the behavior under perturbations 
of the deterministic BS model, introduced in |pj| , |I5| ]. In Sec. M, we study the parallel version of 
it (PBS model), introduced not long ago by Sornette and Dornic [pOj] . We study two versions of 
it, namely the one with random updating and the one with logistic updating. We also present 
here a discussion of the differences and similarities between the BS and the PBS models. In 



Sec. [VI] we discuss the influence of changing the perturbation. Conclusions can be found in Sec. 



VII 



II. GENERAL FORMALISM: DAMAGE-SPREADING IN THE RING 

To illustrate our formalism, let us start by considering a lattice of N sites on a one- 
dimensional ring R\. To each site j we assign a random number fj, extracted from a uniform 



distribution between and 1. We then consider a "replica" R2, in which we introduce a pertur- 
bation by exchanging the positions of the values of fk ± and fk 2 - We define as active the sites k± 
in Ri and (the randomly chosen site) &2 in R2, namely those sites at which we have the same 
value of /. It is clear that, from a statistical point of view, both system R\ and replica R2 are 
described by the same distribution function. The prescription just introduced corresponds, in 
a suitably defined phase space, to a small difference in the initial conditions between R\ and 
i?2 (see Eq. (||) below). Moreover, this procedure of finding an active site and exchanging its 
position with another site taken at random along the lattice, corresponds exactly to the one 



proposed in [14| for the Bak-Sneppen model (other definitions for the initial perturbation are 



considered elsewhere jnj). The dynamics on the ring(s) is defined as follows. At each time-step, 
an integer random number x between 1 and iV is chosen. Bearing in mind that our rings have 
periodic boundary conditions, the positions of the new active sites at time t + 1 is then given 

by 

fci(£ + l) = k x (t)+x (4) 

k 2 {t + 1) = k 2 {t) + x (5) 

on the rings R± and R2 respectively. On these active sites, the value of the variables / is 
changed, assigning to both of them the same random number (this corresponds to the choice 
of the same sequence of random numbers in fl4| or to the same thermal noise in usual damage 
spreading calculations |21~1 ). 

As both system R\ and replica R2 evolve, we compute the Hamming distance, namely 

D(t) = ^p i \f}~f!\- (6) 

Since this quantity has strong fluctuations, we will consider the average (D(t)), over realizations. 
In particular, at t = 1, the average (initial) distance (D(l)) can be obtained from Eq. (^), 

(£>(!)> = |f df l df 2 r ]l U X )ri2{f)\f 1 -f\ (7) 

where rji is the distribution function (at t = 1) for the variable f l e Ri. In this toy model, 
both distributions 77^/) are the same, namely a uniform distribution in /' e R^. A simple 
computation yields 



WD) = 4- (8) 

Applying a similar procedure, one can verify that for 1 ^t ^ N the distance grows linearly. 
Indeed, let us define a(t) as the averaged number of different sites covered in one copy of the 
system at time l^t^N. Then, at time t only a(t) sites have been changed and these are the 
only ones that contribute to distance. From this consideration it follows that 

(D(t)) = (D(l))a(t), (9) 

where the fact that both replica contribute to the distance on the same footing is taken into 
account in Eq. (|8]). In the case of the ring, if N ^> 1 and 1 ^t ^ N the system will always 
choose a new site at each time-step, and therefore a(t) ~ t (note that, in the ID case for this 
dynamics, this is the fastest possible growth of (D(t))). This behavior stops at times t ~ N 
where a crossover to a saturation regime appears. Clearly, after t oc N time steps each site of 
the lattice has been covered at least once. For t ^> N, almost all the lattice sites have been 
covered and the two strings are made of the same random numbers shifted by k2 — k\. Thus, 
the distance reaches a plateau, independent on the size N of the system, given by 

w- 00)) = f'dfdfp^np^nif-fi, (10) 

Jo 
where pi is the normalized distribution function (at t = 00) for the variable f G -Rj. In Eq. 
(|10|), for the particular case of the ring the distributions in the integral are given by pi = rji. 
Applying Eq. (|l(J) to R\ and R% we finally obtain 

\wy(D{t)) = \. (11) 

Note that the same result can be obtained from Eqs. (HH) once o = N is inserted. To have an 
initial distance independent of the lattice size, we consider the ratio (D(t))/(D(1)). For this 
ratio, however, the value of the plateau depends linearly on N, i.e. 

(-D(oo)) N 



(DO)) 



(12) 



In Fig. |I] we show the evolution in time of the ratio LL^ , averaged over many realizations, 
for different lattice sizes. The plateau reached for t — >• 00 depends on iV according to Eq. flT2]). 
The exponent a = 1 obtained for this case from Eq. (^) can also be numerically obtained with 



great accuracy. From a physical point of view, this power law behavior originates in the ability 
of the system to cover the lattice. As we have seen, the activity can jump anywhere on the 
lattice with probability 1/N. Thus the number of sites j with the same fj decreases linearly 
with time and (D(t)) increases linearly with time. As a consequence, the time r needed to 
reach the plateau scales with the lattice size as r ~ N. 

Keeping in mind our goal of modeling the behavior of the BS model, let us now consider 
the case of Levy- Walk type activity jumps along the lattice. More precisely, the length x of 
any jump is extracted from a power law distribution function, namely 

P(x) = (P-l)x-P, (13) 

where the minimum jump is x = 1 and the jump can be to the left or to the right. 

As before, the position of the new active site is obtained by jumping x sites from the present 
one, i.e. the position of the new active site will be given by Eqs. (fUD, where now each copy of 
the system has its own x. Thus, the values of x are un-correlated between the two copies of 
the system. The new values of / assigned to the active sites are the same. This choice results 
in a different behavior at the saturation regime. 

In the Random Walk limit, (3 3> 1 in Eq. (0), the distance Eq. (|7|), can be easily computed 
by considering Eq. @ together with the fact that a ~ t 1 ^ 2 . This calculation yields 

(D(t)} ~ t 1 / 2 . (14) 

In the general case /3 > 1, there is still a power law growth of the distance @ for intermediate 
times 1 ^t ^ N. The exponent in Eq. ([|) can be obtained using the fact that the mean-square 
distance a 2 covered by a Levy Walk behaves like 

t 2 (1< p < 2) 

t 4 ~P (2 < p < 3) 



a 2 {t) ~ < 



(15) 
tlogt (/3 = 3) 



t (/3 > 3) 

in the long-time limit. As a consequence, one can also compute the so-called dynamical expo- 
nent z defined through r ~ N z where r is the time needed for the distance (actually the ratio 



, D(1 1 ) to reach the plateau. This time is given by the time needed to cover all the lattice sites, 
if finite size effects are not counted in. Comparison between Eqs. (0,[|[r5]) yields z = -. 

In this simple model we have excluded any kind of correlation between the values of x 
extracted from Eq. ([H|) and between the two replicas at t = 1. Indeed, the dynamics is given 
by a generalized random walk and therefore the power-law behavior of the growth is not related 
to the statistical properties of the system. This is in fact the idea behind our toy model: We 
use it as a "black box", not knowing what happens inside, we are only able to observe a Levy 
Walk behavior of the activity. Our model is, by conception, a trivial system that has as only 
purpose that of showing what are the consequences, in the context of damage-spreading, of a 
power-law behavior of the activity like the one observed in the Bak-Sneppen model. 

As we shall see in the next section, the non-trivial properties of the self organized critical 
systems are hidden in the value of the exponent a in Eq. (|2|). Furthermore, the averaging 
procedure leading from Eq. (§) to Eq. @ plays a very important role in these non-trivial 
systems. 

Before moving onto the analysis of the BS model, let us discuss in more detail which terms 
are contributing to the computation of (D(t)) via Eq. @. In the ring, we have defined (D(t)) 
by considering the behavior of cr(t), which is a physical quantity related only to the behavior 
of the activity in one single system. In general, considering that the two replica might be 
correlated, we need to update Eq. @ to 

(£>(*)) = (D(l))n cov (t), (16) 

where n cov (t) is the average number of different sites covered by both system and copy. More 
precisely, suppose that at time t, the activity has covered o\ and a 2 different sites in i? x and 
i? 2 respectively. Then, the function n cov (t) is given by 

n CO v(t) = (fli + er 2 - 0i >2 ) . (17) 

where 01.2 represents the number of sites covered in both systems (i.e. the covering overlap 
between system and copy). In the case of the ring, for large N and t <C N, the overlap on the 
rhs of Eq. ([17]) is empty (0^2 — 0) and Eq. (|16D reduces to Eq. (§). In the case of the Bak- 
Sneppen model instead, this intersection cannot be empty even in the thermodynamic limit. 



As a consequence, the exponents predicted from Eq. (p~5|) have to be considered as an upper 
bound for those observables in systems with non-trivial correlations. 

III. THE BAK-SNEPPEN MODEL 

As mentioned above, in its simplest version the BS model describes an ecosystem as a 
collection of N species on a one dimensional lattice. To each species corresponds a fitness 
described by a number / between and 1. For simplicity, one considers periodic boundary 
conditions. The initial state of the system is defined by assigning to each site j a random 
fitness fj chosen from a uniform distribution. The dynamics proceeds in three basic steps: 

1. Find the site with the absolute minimum fitness on the lattice (the active site) and its 
two nearest neighbors. 

2. Update the values of their fitnesses by assigning to them new random numbers from a 
uniform distribution. 

3. Return to step 1. 

After an initial transient that will be of no interest to us here, a non-trivial critical state is 
reached. This critical state, characterized by its statistical properties, can be understood as the 
fluctuating balance between two competing "forces" . Indeed, while the random assignation of 
the values, together with the coupling, acts as an entropic disorder, the choice of the minimum 
acts as an ordering force. As a result of this competition, at the stationary state the majority of 
the fj have values above a certain threshold f c = 0.66702(1) p|. In other words, the distribution 
function of the f/s can be asymptotically approximated by 

Vi(f) = jTT/T ©(/ - /c) , ( 18 ) 

where 0(/) is the Heaviside function. Only a few will be below f c , namely those belonging to 
the running avalanche (see P,P2J for a detailed discussion). Proceeding by analogy with the 



previous cases, once the system is at the critical state, we produce two identical copies Bi and 
_E?2 and find the minimum (the active site). Then, in B2 we swap the value of the minimum 
fitness with the fitness of some other site chosen at random (note that if N is big enough, the 



fitness in the other site will certainly be above threshold). After that, the evolution of the 
Hamming distance given by Eq. (||) is studied. In the evolution of both system and copy the 
same random numbers are used. Here, the length of the jumps in the position of the active 
site follows a power-law distribution given by Eq. ( p~3| ) with (3 ~ 3.23 ||. At variance with the 
case of the ring discussed above, we cannot expect the behavior shown in Eqs. (|9|) and (|T5D. 
Indeed, in the BS model the jumps posses strong spatial correlations, with large probability 
of returning to sites already covered in previous time-steps. As a consequence, the behavior 
of the number of different sites covered in one single system cannot be given by Eq. fll5|) but 



leads to a(t) ~ P* for t > 1 @, with /j, = 0.4114 ± 0.0020 [||]. Moreover, as we consider the 
two copies B\ and B 2 we immediately realize that the two systems are also strongly correlated 
and consequently one obtains an even smaller exponent leading to 



(£>(*)) ~ t a=0 - S2 . (19) 

As we mentioned at the end of Sec. |I|, the behavior of Eq. (|T9| ) can be understood in the 
framework of Eq. flip]). In fact, the decrease in the value of a is given by the appearance of 
a 1>2 ^ in Eq. (|T§. 

Indeed, the fact that we are using the same sequence of random numbers implies that, if 
the system is big enough the absolute minimum in one system is the absolute minimum also 
in the copy. Therefore, if the absolute minimum is among those sites which have not yet been 
covered by the activity, the three terms in the rhs of Eq. ([17D will have the same behavior. 
If the minimum is instead one of the newer values put in the system after perturbation, its 
position on the lattice may be different in the two replica but the three functions in the rhs of 
Eq. ([17]) grow slowly or even do not grow at all. This observation is confirmed by the irregular 
behavior of D(t) in just one single realization. In fact, it is the averaging procedure the one 
that produces finally a smooth curve. It should be noted that, as it is clear from its definition, 
the behavior of the intersection is strongly correlated to the behavior of the other two sets and 
therefore the average in the rhs cannot be split into the sum of the independent averages. As 
a consequence we should expect a smaller exponent with respect to the one obtained for cr(t). 

The initial distance can be computed using Eq. @. To do that we need the distribution 
function of the value of the minimum. Extensive numerical simulations indicate that this 



distribution can be approximated by 

»»(/) = ( 3 -^) 9 ^-/)' ( 20 ) 

where the threshold has been put equal to 2/3. Inserting ( p0|) and fll8l ) in Eq. (||) we obtain 
(£)(!)) ~ 5^. Since (D(t — > oo)) takes into account all sites on the same footing, this saturation 



value can be obtained from Eq. fllOD with pi = P2 = f]\- The distribution 771 comes from Eq. 
(|l8l), and the saturation value is (D(t —* 00)) ~ |. Therefore, as in the case of the ring, the 
saturation value does not depend on the size of the system while the initial distance does. Thus, 
the normalized distance reaches a plateau that must scale with N, as our numerical simulations 
show (see Fig. |2|). 

Coming back to the dynamical exponent z, we find that it still follows the prediction z = 
1/a, as in the case of the ring. For the BS model one obtains, following the above described 
prescription, z ~ 3.12, instead of z ~ 1.6 determined in ||14|| . This value z ~ 3.12 coincides 
reasonably well with the one obtained from the collapse plot (Fig. |3|). The reason for the 
discrepancy between our results and those presented in |TJ| can be traced back to the effects 



of time-rescaling on the (normalized) growth function. Indeed, let us assume we use a different 
time-scale, and consider the case in which we make a measure of D{t) every v time-steps 
(instead of every time-step), where v is distributed according to a certain function P{y). The 
rescaled distance D{t) will be given by 

D(t) = fdvP(v)D((t-l)(v) +u) , (21) 

where (v) = J dvvP(v) is the average number of time-steps between two consecutive measure- 
ments. The choice made in |TJ]] corresponds to P(v) = 8{y — N). It is easy to see that, in this 
case, the growth exponent for D(t) is still a, but the measured dynamical exponent is given by 
1/a — 1. In principle, one could imagine more complicated distributions P(y) for the measur- 
ing time. In particular, if P(v) did not have a finite first moment (as would be the case, for 
instance, if P{v) corresponded to the avalanche distribution) Eq. (|2l|) would yield D(2) ~ N, 
i.e. the rescaled distance would saturate almost immediately |24[ . 



In order to study the size dependence of a we have performed a set of measurements for 
different system sizes and then extrapolated the value of a to infinite size. We get an asymptotic 



value a asym = 0.40(1) (see Fig. |). 

At this point, it is worth discussing what happens in higher dimensions. The high- 
dimensional Bak-Sneppen model has been extensively studied in [^j , where the behavior of the 
exponent /i for the growth of the quantity a(t) has been computed until the mean- field regime 
/i = 1 has been reached. In the framework of the damage-spreading, taking into account the 
correlations as discussed above, one expects the exponent a to follow a similar pattern, reaching 
the value a = 1 in the mean-field case. These mean-field results can also be obtained in the 
random- neighbors case ^6f . However, from the point of view of damage spreading, the random 



nearest neighbor case presents a complication. There is an ambiguity in the choice of the neigh- 
bors. Indeed, their absolute positions on the lattice should be either the same in the two copies 
of the system or taken at random in an uncorrelated fashion. In both cases, each one of the 
two copies will behave normally, but the behavior of the distance will be completely different. 
Indeed, in the first case the distance between the two systems will never grow, while in the 
latter case the behavior of the distance resembles that of the ring with uniformly distributed 
jumps. 

IV. BS MODEL WITHOUT NOISE 



Through the use of different maps (chaotic as well as non-chaotic), it was shown flq , |T9|j that 
the random updating is not a necessary requirement to have SOC Moreover, as long as the 
updating rule is chaotic the system does not change the universality class, i.e. all the exponents 
are the same as in the case of random updating. This means that the system is able to self- 
organize at a higher level: It takes into account the temporal correlation (or the average time 
spent in every site) by increasing the threshold, so as to have the same statistical properties 
|I~8| , [T9"|| . As a consequence, all equations and relations derived for the original BS model are 



still valid for all the cases with chaotic updating. The stationary distribution of the fitnesses, 
on the other hand, follows a different pattern. Indeed, the position of the threshold as well as 
the exact shape of the stationary distribution depends on the actual form of the updating rule. 
To calculate the value of (-D(l)) one needs to specify the nature of the perturbation. Using 
the same methods as in the random updating case, i.e. to swap the position of the minimum 



in B2 with any site taken at random, at t = 1, the average (initial) distance (D(l)) can be 
obtained from Eq. (0), where r\i is the distribution function (at t — 1) for the variable /* G B ri . 
When one considers the chaotic updating, Eqs. (|T|,|10D are still valid. The only difference 
is that one needs to find the stationary distributions that correspond to the actual map. For 
instance, for the tent map as well as for the Bernoulli map the distribution coincides with that 
of the random updating (except for the value of the threshold in the Bernoulli case) fll8"Ji~9 ]. As 
an example, we substitute the random updating with the logistic map 

/i(t+l) = 6/i(*)(l -/<(*)), (22) 

where % runs over the minimum and its nearest neighbors and b is a parameter (that we set to 
the value 4 which is at the threshold of the chaotic phase). Inserting the estimated distributions 
obtained in [[H|[n| one obtains as first approximation D asym ~ 0.087 and the same exponent a 
as in the random updating. 

In Fig. [5] one can see the evolution of the distance for the case of a Bak-Sneppen model 
with random and logistic updating, using the perturbation defined in ]TJJ (flipping of the 
minimum fitness with another fitness chosen at random). Simulations with other chaotic maps 
give the same exponents too: This reflects the observed fact that chaotic maps do not change 
the universality class of the model J18|,|l9| . In the inset we show the scaling of ND asym (N) 
versus N for both random and logistic updating. A linear fit gives as estimation of the plateau 
Dasym = 0.112(5) for random updating and D asym = 0.088(5) for logistic updating, in good 
agreement with our analytic estimate. 

V. PARALLEL BS MODEL 



The (i-dimensional parallel version of the BS model (PBS) |2C|| , where parallel updating has 
been introduced, has been found to have an exact mapping onto 2d Directed Percolation. In 
fact, the avalanche time distribution in the PBS model is equivalent to the cluster distribution 
in DP and the threshold for PBS is equivalent to the critical probability in DP. 

The dynamics in the PBS is different from the one in the (usual) Bak-Sneppen model due 
to the introduction of parallel updating: 



1. Find the site with the absolute minimum fitness f m i n on the lattice (the active site) and 
its two nearest neighbors. 

2. Update the values of their fitness by assigning to them new random numbers from a 
uniform distribution. 

3. Search for all sites on the lattice with fitness / < f m in and update them together with 
their nearest neighbors. Repeat the search until there are no sites left with / < f m i n . 

4. Return to step 1. 

The difference between the PBS model and the original BS model lies in step 3. In the 
extremal version, once the minimum and its nearest neighbors are updated one looks for the 
new minimum, and consequently the number of updated sites per time step, U t , is constant 
(Ut = 3). In the parallel version, instead, this number will follow a complex temporal evolution 
during an avalanche with, in general, Ut > 3. In fact, the distribution of the number of updated 
sites per time step (inside an avalanche) shows a nearly flat distribution, with a upper cutoff, 
whose value is comparable with the system size |17| . Due to this saturation, one observes that 
D(t) grows in time faster than in the BS case and that finite size effects are also much stronger. 

Since the disorder is stronger in the parallel version (Ut > 3), one expects the equilibrium 
point to be displaced towards the completely disordered value. In fact, for the PBS model, 
f c ~ 0.5371(1) p0[ . This results can be obtained both numerically and analytically by mapping 
the model onto directed percolation. 

To study the behavior under perturbations of the PBS model, we follow the same procedure 
as for the extremal case. We produce two identical copies B\ and B 2 of the system of size 
N in the critical state, and find the minimum (the active site). Then, we introduce a slight 
perturbation in B 2 and follow the evolution of the Hamming distance Eq. (^) in time. Since 
this quantity has strong fluctuations, we consider the average (D(t)) over different realizations 
of the initial values of the fj. In particular all the simulations presented here, are the result 
of averaging over 10 2 realizations. The simulations are performed with both random and 
deterministic (logistic) updating rules. Let us begin by discussing the results obtained for 
random updating. As mentioned above, D(t) may depend on the internal correlation of the 



system and on the correlations between the two copies. In the ID BS model, due to the choice 
of unit of time, which allows only a number 0(1) of sites to be updated, the growth rate must 
give an exponent a < 1 and stop at a certain time r ~ N z , at which a crossover to a saturation 
regime appears. 

In the PBS case, too, D(t) reaches, after an initial power law growth (as in the BS case), 
a well defined plateau (see Fig. |6|). However, due to the faster parallel dynamics, the value of 
the exponent a is found to be larger than the one for the the extremal non-parallel case, as 
shown in Table |. Moreover, there seem not to be dependence of the exponent a on the system 
size N. Notice that our exponent differs from the one obtained in |27|] for the Domany-Kinzel 



model in the context of DP, onto which the PBS can be mapped. This discrepancy is due to 
the choice of time-scale for the measure of the distance. Indeed, in the Domany-Kinzel model 
case only one active (occupied) site per time step can be updated, together with its neighbor, 
thus resulting in a dilated time scale with respect to the study presented here. Thus, in order 
to compare the two models one has to establish the relationship between the two time-scales. 
This is not easy to realize for the Hamming distance. In fact, according to this interpretation, 
equal times for the two copies on the Monte-Carlo parallel time-scale are not equal times on 
the DP-like time-scale. To circumvent this problem, we realized a set of PBS simulations (with 
system size N = 2000) for a single copy, computing at every Monte-Carlo (parallel) time step 
St — 1 the number n act (t) of sites below threshold (which is itself time-dependent). This defines 
the temporal increment for the DP-like time scale Stop = St-n act (t) = n act (t). The effective DP 
time at the MC step t is thus connected to the effective time at MC step t — 1 by the relation: 

t D p(t)=t DP (t-l)+n act (t). (23) 

Then, we mediated over different realizations of the dynamics, obtaining the scaling of the 
effective DP time with the MC time of the simulation. The result, shown in Fig. |j], is that 
top ~ t*> with C, = 1.41(2). From this, if we assume that top is the equivalent of the DP time- 
scale for PBS, we can combine the scaling law for D(t) and that for top to get the effective 
scaling exponent a* for the Hamming distance with respect to the DP time-scale top'- 

= a (L47 = a33(1) _ / 24 x 

C 1.41 V ; K ' 

This value is quite near to the DP exponent aop = 0.314(1) 



At any given time step t during an avalanche, the average growth of D(t) is connected to 
the mean number of sites a(t) covered by the activity in each system according to Eq. (|9]). In 
Fig. [8], a(t) exhibits a power law behavior, a(t) ~ t M . The values of the scaling exponent jjl 
at different sizes N are always bigger than the corresponding values of the exponents for the 
Hamming distance (Tab. []]). This is expected since the correlations between the two copies in 
D(t), if present, can only decrease the value of a with respect to /i [jl5f . The asymptotic value 
of the Hamming distance shows quite a strong and persistent dependence on the system size N 
and converges to an asymptotic value only logarithmically in N (see inset of Fig. |6|). Then, in 
order to get the real value of the plateau, one has to go to the thermodynamic limit N = oo, 
once the plateau has been reached. This is realized by a logarithmic extrapolation of the 
data for different sizes. The value of this plateau can be obtained in terms of the asymptotic 
stationary fitness distribution p(f) from Eq. ( |T0D where Pi{f) is the normalized distribution 
function (at t = oo) for the variables /* 6 -Hj, at a given system size N. The dependence of the 
plateau D asym (N) on N must then be related to the shape of the stationary fitness distribution. 
The fitness distribution has, in fact, a flat tail below the threshold f c , which disappears only 
logarithmically in N as the system size is increased. In the limit N — > oo the distribution 
Pi(f) is given by Eq. flTSP with f c = 0.5371. By substituting the distribution pi thus obtained 
into Eq. ( fLOf) we get D asym (N = oo) = D asym = 0.1543 as estimation of the plateau in the 
thermodynamic limit. Turning back our attention to Fig. ^|, one can see that this result is 
consistent with the extrapolated numerical value D asym ~ 0.14(2). 

We have also performed a similar analysis for the PBS model with a deterministic updating, 
in which the new fitnesses are obtained by iterating the logistic map, namely 

fi(t + l) = bfi(t)(l-Mt)), (25) 



where fi(t) is the fitness of site i at time t, and b is a parameter set to the value 4 f!l|[T9"|j . We 
observe the same qualitative behavior for all the quantities studied in the random updating case. 
The fitness distribution is however different since it is influenced by the invariant measure of the 
map, as pointed out in [|l8i|l9l . In the present case, the fitness distribution is strongly picked 
around / = and around / = 1, for finite sizes N. The distribution is of course not symmetric 



and in the large N limit, all the /,• are above a threshold f c = 0.55(1) p8| . The value of the 



plateau D asym converges, in the limit N — > oo, to D asym ~ 0.15(2). By inserting the fitness 
distributions thus obtained in Eq. (|T0|), we obtain an analytic estimation, D asym = 0.1556, 
that is indeed very close to the random updating analytical value, and in agreement with our 
numerical results. This is reasonable, since the threshold of the parallel BS with deterministic 
rule is very near the threshold of the random updating case. Although we performed, for the 
computation of the plateau, simulations up to system size N = 8000, the need of a logarithmic 
extrapolation towards N = oo prevents us from obtaining a precise numerical estimation of the 
plateau. The values of the exponents a and \x show no substantial differences with respect to 



the random updating case, as is the case for the extremal BS model ||16| , thus confirming the 
robustness of a with respect to different updating rules. 

VI. ON THE INITIAL PERTURBATIONS 

One point remains, however, that needs to be studied. The definition of the initial pertur- 



bation in the replica in |H] is too restrictive. In fact, by considering as initial perturbation 

fi = fi + egi , (26) 

where e is a small positive number, we can take several choices for g t without altering the 
exponent a. We considered four different implementations of go 

(a) gi = ip(t) where ip(t) is a random number between and 1. 

(b) gi = ip(t) where ip(t) is a random number between -1/2 and 1/2. 

(c) gi = ip(t) for i corresponding to the minimum and zero otherwise. 

(d) gi = constant = 1 for all sites. 

This kind of perturbation allows us to tune the initial mean distance D(t), to any arbitrarily 
small value depending on e in Eq. Q2SD (contrary on the flipping introduced in |TJ| , which gives 



a fixed, size dependent, mean initial distance). In case (d), for example, (which is the case 
we will use in the analysis below) the initial distance is D(l) = e, independent on N, and the 
plateau will be independent on iV too. This characteristic of the global perturbation is useful 



to explore some properties of the model with deterministic updating and is more in line with 
the standard techniques of damage spreading problems. 

We have performed the analysis with the prescription ( p6|) also in the case in which the model 
has a deterministic microscopic rule, like an updating with the logistic map. In this case, the 
numbers in each replica will not be the same since the map will be applied on different numbers. 
In Fig. P we show the behavior of D(t) for random and logistic updating, where we used 
perturbation (|26| ) with the implementation (d). The exponent does not change: In both cases a 
value of a around the value 0.32 found in |L4| is obtained. The plateau of ND(t) (we computed 
ND(t) to reduce statistical fluctuations) for the case of global perturbation, scales linearly 
with the system size N, as shown in the inset of Fig. [7|, with a slope a = D asym = 0.114(5) 
and a = D asym = 0.089(5) respectively for random updating and logistic updating. The good 
agreement with the rough analytic estimate obtained above and with numerical results for the 
flipping shows that the value of the plateau, too, does not depend on the kind of perturbation 
applied. Consequently, the distance D{t) has a plateau D asym independent on the size, as it 
happens with the flipping. These results do not change if we use the implementations (a), (6), (c) 
of the perturbation (^). 

The case of the logistic map with perturbation (12(1) is particularly interesting from a dif- 
ferent point of view. Since the map is chaotic, one could expect that on small scale distances 
(-D(i) <C 1/-/V) the chaoticity of the maps dominates and the distance of the two copies grows 
exponentially. For bigger distances instead, the dynamics is dominated by the damage spread- 
ing (and thus the critical properties of the BS model) and the distance grows as a power law. 
This is exactly what we find numerically. In Fig. [H] we show numerical simulations of the BS 
model with deterministic updating with the logistic map (with parameter value 4) and random 
perturbation, for different system sizes and initial distances. The distance D(t) has an initial 
exponentially growing phase with a Lyapunov exponent A = 0.49(1). This Lyapunov exponent 
is smaller than the Lyapunov exponent of a single logistic map with 6 = 4, which actually is 
log2 = 0.693147.... Indeed, the microscopic dynamical rule of the system changes the values of 
the minimum fitness and its nearest neighbors, and the Lyapunov exponent we measure arises 
from the interaction between the three maps applied to the three different numbers. 



VII. CONCLUSIONS 

We have shown that the power-law behavior of the distance Eq. (Q) has its origins in the 
behavior of the mean-squared distance covered by the activity. The first consequence is that 
a < 1. The second is that the internal correlations of the jumps, governed by Eq. (|T3|), and the 
strong correlations between the two copies, can severely slow down the growth of (D(t)), thus 
leading to smaller exponents with respect to the ones predicted in Eq. (|T5|). Since an analytic 
derivation of the exponents characterizing the critical properties of the BS model is yet to be 
found, we had to base our work on numerical results. Nevertheless, by rewriting Eq. (£J) into 
the more appropriate form given by Eqs. (|l6l - |l7|) , we have been able to study the mechanisms 
leading to it. In the framework of these equations the occurrence of the plateau can be also 
easily understood and we have computed its value. 

Our theory assumes that neither in the BS case nor in the case of the ring will the distribution 
of the variables / change during the measurement of (D(t)) (we exclude the transient). 

Moreover, we have studied the characteristics of damage spreading in a parallel version of 
the Bak-Sneppen model, and compared them to those of the extremal version. In general, the 
PBS model exhibits a faster evolution towards a stationary state, which is strongly dependent 
on the system size N. 

Finally we have confirmed the picture presented in |nj according to which the application 
of deterministic (chaotic) maps does not change the universality class of the model. Further- 
more, we have shown here that the analysis attempted in |L4| and |L5| is independent of the 
microscopic details of the system, and of the prescription of the initial perturbation. This of 
particular importance if one looks at these results within the framework of standard damage- 
spreading. The possibility to apply this analysis to such different versions of the Bak-Sneppen 
model, together with the good agreement with the results for damage- spreading in Directed 
Percolation, gives us confidence that these techniques can be extended to other Self- Organized 
Critical models. 
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TABLES 



N 


500 


750 


2000 


4000 


8000 


a ran. up. 


0.48(1) 


0.48(1) 


0.47(1) 


0.49(1) 


0.48(1) 


a log. up. 


0.47(1) 


0.47(1) 


0.49(1) 


0.48(1) 


0.48(1) 


fj, ran. up. 


0.63(1) 


0.63(1) 


0.65(1) 


0.65(1) 


0.65(1) 


H log. up. 


0.63(1) 


0.64(1) 


0.63(1) 


0.65(1) 


0.65(1) 



TABLE I. Values of the exponents a and \x for different sizes and updating rules. 
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FIG. 1. Distance D(t) between the two replicas for three different sizes N for the uncorrelated 
ring model. 



][)[K) 

N 



FIG. 2. Scaling of the long time plateau D* = /^q-a as a function of the number of sites N for 
the BS model. The best fit yields an exponent of 1.03(4). 




FIG. 3. Collapse plot of the data for the evolution of the distance, in the BS case. The best fit 
yields a = 0.32(3) and z = 3.0(2). 
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FIG. 4. Extrapolation of a for the BS model with random updating. 
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FIG. 5. log 10 -log 10 plot of D(t) for the BS model with flipping, with random updating and de- 
terministic updating with logistic map, for different system sizes. The value of the plateau actually 
depends on the updating rule, but it does not depend on the system size. Inset: A plot of the plateau 
of ND(t) versus N fits very well with a linear scaling, showing that the plateau of the Hamming 
distance D{t) is size independent. 
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FIG. 6. The logio — logio plot of the Hamming distance D(t) in the parallel BS model, versus 
the time t. Inset: logarithmic extrapolation of the plateau as a function of N, for both random and 
logistic updating rules. The infinite size value D asym obtained is in agreement with our analytical 
estimate. 
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FIG. 7. Scaling of the effective DP-like time top with the Monte-Carlo time t, on a logio ~ l°9w 
scale. We get a power law behavior with exponent £ = 1.41(2). 
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FIG. 8. logio ~ l°9w plot of the number a(t) of sites covered during the evolution of a single copy 
of the system, versus the time t, for random updating rule. 
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FIG. 9. logxQ-logxo plot of D{t) for the BS model with the perturbation (pq), implementation (d), 
with random updating and deterministic updating with logistic map, for different system sizes. The 
value of the plateau actually depends on the updating rule, as in the flipping case. Inset: A plot of 
the plateau of ND(t) versus iV gives the same findings of the case of flipping. 
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FIG. 10. Zmear-log 10 plot of D{t) for the BS model with perturbation (f2q), implementation (d), 
and updating with a logistic map with parameter 6 = 4, for different system sizes and mean initial 
distances e. 



